{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example: Vertical instability (growth rate) calculations\n", "\n", "---\n", "\n", "This notebook will demonstrate various different ways to calculate the **growth rates** (or **timescales**) associated with **vertically unstable modes** of an equilibrium. We will explore ways to plot these different quantities, carry out eigenmode reduction, and carry out more advanced analysis. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "### The deformable growth rate model\n", "\n", "Assume we have solved the Grad-Shafranov equation and have obtained an equilibrium. From this we can extract the plasma current density discretised over the spatial domain (whose points are only non-zero inside the LCFS, so we mask points outside the limiter to save memory):\n", "\n", "\\begin{equation*}\n", " \\vec{I}_y := \\vec{I}_y^{GS} (\\vec{I}_m, I_p, \\vec{\\theta}).\n", " \\tag{1}\n", "\\end{equation*}\n", "\n", "This is a function of the: \n", "- currents in the metal conductors $\\vec{I}_m$ (i.e. the active coils and passive structures).\n", "- plasma current $I_p$.\n", "- parameters $\\vec{\\theta}$ used to define the plasma current density profiles (i.e. $p'(\\psi_n)$ and $FF'(\\psi_n)$). \n", "\n", "In order to estimate the instability timescales/growth rates of this plasma, we will need to define (and then linearise) the coupled **circuit** and **plasma** equations. \n", "\n", "**Circuit equations**: \n", "\n", "These equations govern the flow of current in the active coils and passive structures (while being coupled to the currents in the plasma):\n", "\n", "\\begin{equation}\n", " M_{m,m} \\vec{\\dot{I}}_m + M_{m,y} \\vec{\\dot{I}}_y + R_{m,m} \\vec{I}_m = \\vec{V}_m\n", " \\tag{2}\n", "\\end{equation}\n", "\n", "- $M_{m,m}$ is the symmetric matrix of mutual inductances between all of the metals. \n", "- $M_{m,y}$ is the matrix of mutual inductances between all of the metals and the discretised plasma current density. \n", "- $R_{m,m}$ is the diagonal matrix of resistances in the metals.\n", "- $V_m$ is the vector of voltages applied to the metals (note: these are non-zero for the active coils only).\n", "- Note that the dots indicate time derivatives. \n", "\n", "**Plasma equation**: \n", "\n", "This equation governs the flow of current in the plasma (while being coupled to the currents in the metals):\n", "\n", "\\begin{equation}\n", " \\frac{\\vec{I}_y^T}{I_p} \\left( M_{y,y} \\vec{\\dot{I}}_y + M_{y,m} \\vec{\\dot{I}}_m + R_{y,y} \\vec{I}_y \\right) = \\vec{0}\n", " \\tag{3}\n", "\\end{equation}\n", "\n", "- $M_{y,m} = M_{m,y}^{T}$ (from above).\n", "- $R_{y,y} = 2 \\pi \\sigma_p \\hat{R} /dA$ is the diagonal matrix of resistances for each plasma element.\n", " - $\\sigma_p$ is the assigned plasma resistivity. \n", " - $\\hat{R}$ is a matrix storing the radial position of each plasma element.\n", " - $dA$ is the area of each plasma element.\n", "\n", "\n", "**Linearisation of discretised plasma current density**:\n", "\n", "To simplify these equations, we linearise equation (1) around the current GS equilibrium such that:\n", "\n", "\\begin{equation}\n", " \\dot{\\vec{I}_y} \\approx \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS} \\dot{\\vec{I}}_m + \\frac{\\partial \\vec{I}_y}{\\partial I_p}\\Bigg|_{GS} \\dot{I}_p + \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}}.\n", " \\tag{4}\n", "\\end{equation}\n", "\n", "The partial derivative terms (i.e. the Jacobians) are calculated in FreeGSNKE using finite differences. The use of these Jacobians is what makes this a **deformable** growth rate model. \n", "\n", "We can plug this equation into both (2) and (3) and re-arrange:\n", "\n", "\\begin{align}\n", " \\left( M_{m,m} + M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS} \\right) \\vec{\\dot{I}}_m + M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial I_p}\\Bigg|_{GS} \\dot{I_p} + R_{m,m} \\vec{I}_m + M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}} = \\vec{V}_m \\tag{5} \\\\\n", " \\frac{\\vec{I}_y}{I_p} \\left[ \\left( M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS} + M_{y,m} \\right) \\dot{\\vec{I}}_m + M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial I_p}\\Bigg|_{GS} \\dot{I_p} + R_{y,y} \\vec{I}_y + M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}} \\right] = \\vec{0}. \\tag{6}\n", "\\end{align}\n", "\n", "Next, we need to write this in state matrix form. For this, we need to remove the dependence of equation (6) on the denominator $I_p$ by introducing two terms:\n", "- $\\hat{\\vec{I}}_y = \\vec{I}_y / I_p$ is the normalised discretised plasma current density over the spatial domain.\n", "- $R_p = \\hat{\\vec{I}}_y^{T} R_{y,y} \\hat{\\vec{I}}_y$ is the \"lumped\" plasma resistance.\n", "\n", "This results in the state matrix system $\\vec{M} \\dot{\\vec{x}} + \\vec{R} \\vec{x} = \\vec{c}$:\n", "\n", "\\begin{equation*}\n", "\\begin{aligned}\n", "\\left[\\begin{array}{cc}\n", " M_{m,m} + M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS} & M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial I_p}\\Bigg|_{GS} \\\\ \n", "\\frac{\\hat{\\vec{I}}_y}{R_p} \\left( M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS} + M_{y,m} \\right) & \\frac{\\hat{\\vec{I}}_y}{R_p} M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial I_p}\\Bigg|_{GS}\n", "\\end{array}\\right]\n", "\\left[\\begin{array}{c}\n", "\\dot{\\vec{I}}_m \\\\ \n", "\\dot{I}_p\n", "\\end{array}\\right] +\n", "\n", "\\left[\\begin{array}{cc}\n", "R_{m,m} & 0 \\\\ \n", "0 & 1\n", "\\end{array}\\right]\n", "\\left[\\begin{array}{c}\n", "\\vec{I}_m \\\\ \n", "I_p\n", "\\end{array}\\right]\n", "\n", "&=\n", "\\left[\\begin{array}{c}\n", "\\vec{V}_m \\\\ \n", "0 \n", "\\end{array}\\right]\n", "-\n", "\\left[\\begin{array}{c}\n", "M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}} \\\\ \n", "\\frac{\\hat{\\vec{I}}_y}{R_p} M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}}\n", "\\end{array}\\right].\n", "\n", "\\end{aligned}\n", "\\end{equation*}\n", "\n", "Re-arranging, we can write this as \n", "\n", "$$ A\\dot{\\vec{x}} = \\vec{x} - \\vec{b}, $$\n", "\n", "where $A = -\\vec{R}^{-1} \\vec{M}$ and $\\vec{b} = \\vec{R}^{-1} \\vec{c}$. \n", "\n", "**Growth rates**:\n", "\n", "To ascertain the growth rates of the system, we need to look at the eigenvalues of the system in the absence of any control voltages and change in profile parameters (i.e. $\\vec{c}=\\vec{0}$). The eigenvalues of $A$, denoted as $\\tau_i$, represent the **timescales** of the system (units in seconds). The corresponding eigenvectors $\\vec{v}_i$ are the **modes** of the system and represent currents in the metals and the total plasma (the last element).\n", "\n", "Therefore:\n", "- $\\tau_i > 0$ corresponds to an **unstable** mode $\\vec{v}_i$.\n", "- $\\tau_i < 0$ corresponds to a **stable** mode $\\vec{v}_i$.\n", "\n", "These timescales describe the evolution of the plasma vertical instability on the resistive timescale of the conducting metals. Its reciprocal, $\\gamma_i = 1/\\tau_i$, gives the **growth (or decay) rate** (units in 1/seconds) of a particular mode (decay rate if $\\gamma_i < 0$). \n", "\n", "These rates tell us how fast a perturbation in the system grows or decays and later on we will visualise these (eigenmode) perturbations by plotting the poloidal flux they produce in the respective metals.\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Basic usage\n", "\n", "Here we'll use the same equilibrium as in previous examples, except we've slightly modifed the coil currents to generate a more elongated (vertically unstable) plasma. These coil currents were found by shifting the lower X-point vertically downwards using technqiues from the virtual circuits example notebook." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "import matplotlib.pyplot as plt\n", "import freegs4e\n", "import numpy as np\n", "\n", "# build machine\n", "from freegsnke import build_machine\n", "tokamak = build_machine.tokamak(\n", " active_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_active_coils.pickle\",\n", " passive_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle\",\n", " limiter_path=f\"../machine_configs/MAST-U/MAST-U_like_limiter.pickle\",\n", " wall_path=f\"../machine_configs/MAST-U/MAST-U_like_wall.pickle\",\n", ")\n", "\n", "# initialise equilibrium object\n", "from freegsnke import equilibrium_update\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak,\n", " Rmin=0.1, Rmax=2.0, # radial range\n", " Zmin=-2.2, Zmax=2.2, # vertical range\n", " nx=65, # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", " # psi=plasma_psi\n", ") \n", "\n", "# initialise profile object\n", "from freegsnke.jtor_update import ConstrainPaxisIp\n", "profiles = ConstrainPaxisIp(\n", " eq=eq,\n", " paxis=8.1e3,\n", " Ip=6.2e5,\n", " fvac=0.5,\n", " alpha_m=1.8,\n", " alpha_n=1.2\n", ")\n", "\n", "# initialise solver\n", "from freegsnke import GSstaticsolver\n", "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) \n", "\n", "# set coil currents (these were found using the virtual circuits example)\n", "current_values = {'Solenoid': 5000.0,\n", " 'PX': 4664.8264042905585,\n", " 'D1': 4817.449859117131,\n", " 'D2': 1234.1293690145924,\n", " 'D3': 1369.133338804566,\n", " 'Dp': -905.4445404924389,\n", " 'D5': 3757.2659244120455,\n", " 'D6': -675.5868544077532,\n", " 'D7': -94.84095712236905,\n", " 'P4': -3536.0809826602153,\n", " 'P5': -4568.000388616427,\n", " 'P6': 0.0005793038603228533,\n", "}\n", "\n", "for key in current_values.keys():\n", " eq.tokamak.set_coil_current(coil_label=key, current_value=current_values[key])\n", "\n", "# carry out forward solve\n", "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=None, \n", " target_relative_tolerance=1e-9)\n", "\n", "# plot the resulting equilbrium\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=60)\n", "ax1.grid(True, which='both')\n", "eq.plot(axis=ax1, show=False)\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solving the linearised system requires inititialising the `nl_solver` object. This requires the:\n", "- equilibrium and profiles. \n", "- plasma resistivity (for the lumped plasma resistance calculation).\n", "- some options regarding mode cut-off frequencies, minimum Jacobian column norms, etc (we will revist these later in the notebook so we can ignore these for now).\n", "\n", "Within this object, the Jacobians are calculated and the matrices are built. Depending on the number of active coils and passive structures present in the system, this can take a while (as a GS solve is required for each). Following this, the matrices of the linearised system are used to find the eigenvalues/vectors.\n", "\n", "Here we will be using **all** available metal modes, however, later on we will show how to exclude those (via additional options to `nl_solver`) with specific frequencies, timescales, coupling, etc to reduce computational runtime (without losing too much information)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import nonlinear_solve\n", "\n", "nonlinear_solver = nonlinear_solve.nl_solver(\n", " eq=eq, \n", " profiles=profiles, \n", " GSStaticSolver=GSStaticSolver, \n", " plasma_resistivity=1e-6, # this defines the lumped plasma resistances\n", " min_dIy_dI=0, # this has been set artificially low in this example\n", " threshold_dIy_dI=1, # this has been set artificially high in this example\n", " max_mode_frequency=1e10 # this has been set artificially high in this example\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the printout we can see the steps being taken.\n", "\n", "Once instatiation has taken place, the \"mode selection criteria\" is chosen based on user inputs - again, we revist this later. This is followed by the \"initial mode selection\" We see that all active coils (12) and all passive structure modes (138) are included in this case.\n", "\n", "Following this, the Jacobian is built. This has dimensions \"number of plasma current density points inside the limiter\" x \"currents included in the analysis\". The currents in the analaysis are the active coils (12), the passives (138), and the total plasma current (1) from the plasma equation (total modes being 151). \n", "\n", "After this, any further mode reduction takes place (in this case none) and the Jacobian is re-sized if required. \n", "\n", "The 'Stability parameters' section shows a number of different metrics used to gauge plasma vertical instability. In 'Deformable plasma metrics' we have the growth rate of the unstable mode and the associated timescale (i.e. 1/(growth_rate)). In the 'Rigid plasma metrics' section, we return a number of commonly used parameters that do not make use of the Jacobian matrix (we will discuss these later too)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Extracting data and visualising\n", "\n", "Next, we will show to access some of the calculated data and how to visualise a few things. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# accessing the number of modes\n", "print(f\"Total number of modes excl. plasma current = {nonlinear_solver.n_metal_modes}\") # total (actives + passives)\n", "print(f\"Total number of active coils = {nonlinear_solver.n_active_coils}\") # actives\n", "print(f\"Total number of passive structures = {nonlinear_solver.n_passive_coils}\") # passives" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# accessing the growth rates (via the timescales)\n", "timescales = nonlinear_solver.linearised_sol.all_timescales # all eigenvalues: timescales\n", "growth_rates = 1/timescales # growth rates are simply 1/timescales\n", "modes = nonlinear_solver.linearised_sol.all_modes # all eigenvectors (columns in same order as e'values)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extracting the unstable mode\n", "mask = (timescales > 0)\n", "idx = np.where(mask)[0][0] # index of unstable mode\n", "unstable_timescales = timescales[mask]\n", "unstable_modes = np.squeeze(modes[:,mask])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can visualise the poloidal flux produced by the unstable mode using the following code. We can also do this for other (stable) modes if required (just change the index `i` in the cell below)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# mode number (choose which one you want to visualise)\n", "i = idx # default is unstable mode\n", "mode_currents = np.real(modes[:,i])\n", "\n", "# the associated instability timescale and growth rate\n", "print(f\"Mode {i} ---> {'stable' if np.real(timescales[i]) < 0 else 'unstable'}\")\n", "print(f\"Growth rate = {np.real(growth_rates[i]):.2e} [1/s]\")\n", "print(f\"Timescale = {np.real(timescales[i]):.2e} [s]\")\n", "\n", "# multiply each metal current (from the eigenvector) with its corresponding Greens matrix and sum\n", "# (don't forget to omit the plasma current mode, i.e. the final element)\n", "flux = np.sum(mode_currents[0:-1, np.newaxis, np.newaxis] * nonlinear_solver.vessel_modes_greens, axis=0)\n", "\n", "# plot\n", "fig, ax = plt.subplots(1, 1, figsize=(5, 8), dpi=60)\n", "ax.grid(True, which='both')\n", "ax.set_aspect('equal')\n", "ax.set_xlim(0.1, 2.15)\n", "ax.set_ylim(-2.25, 2.25)\n", "ax.set_title(f\"Poloidal flux (mode {i})\")\n", "ax.set_xlabel(r'Major radius, $R$ [m]')\n", "ax.set_ylabel(r'Height, $Z$ [m]')\n", "plt.tight_layout()\n", "\n", "eq.tokamak.plot(axis=ax,show=False)\n", "ax.plot(eq.tokamak.wall.R, eq.tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im = ax.contour(eq.R, eq.Z, flux, levels=50) \n", "cbar = plt.colorbar(im, ax=ax, fraction=0.09)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we can extract two different plasma current density $J_{p}(\\psi, R,Z)$ maps: call them $J_{p}^{deform}$ and $J_{p}^{rigid}$.\n", "\n", "The first corresponds to an application of the unstable mode currents to the equilbirium, solving the GS equation, and observing the \"deformable\" movement of the plasma core (it shifts in $(R,Z)$ and also the boundary will displace). This results in a new current density map $J_{p}^{deform}$. \n", "\n", "The second corresponds to taking the $(R,Z)$ shifts in the current centre (from $J_{p}^{deform}$) and shifting the $J_p$ map from the orginal equilibrium. In this case we observe a \"rigid\" displacement of the plasma current density: $J_{p}^{rigid}$. \n", "\n", "While they do not look significantly different in the case presented here, they will when studying more vertically unstable plasmas. Below we plot the boundaries of these two Jtor maps (left) and the difference between the full maps (right). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# data\n", "jtor_maps = nonlinear_solver.deformable_vs_rigid_jtor \n", "\n", "# rate of change of R and Z current centre of plasma wrt to the unstable mode only\n", "dRZd_unstable_mode = nonlinear_solver.dRZd_unstable_mode\n", "print(f\"Rate of change of Rcurrent wrt unstable mode = {dRZd_unstable_mode[0]:.2e} [m].\")\n", "print(f\"Rate of change of Zcurrent wrt unstable mode = {dRZd_unstable_mode[1]:.2e} [m].\")\n", "\n", "diff = np.abs(jtor_maps[0] - jtor_maps[1])\n", "\n", "# plot\n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 8), dpi=60)\n", "\n", "ax1.grid(True, which='both')\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_title(f\"Deformable and rigid boundaries\")\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')\n", "plt.tight_layout()\n", "\n", "eq.tokamak.plot(axis=ax1,show=False)\n", "ax1.plot(eq.tokamak.wall.R, eq.tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "ax1.contour(eq.R, eq.Z, jtor_maps[0], levels=[0], colors='b') \n", "ax1.contour(eq.R, eq.Z, jtor_maps[1], levels=[0], colors='r') \n", "\n", "\n", "ax2.grid(True, which='both')\n", "ax2.set_aspect('equal')\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "ax2.set_title(f\"Abs. diff. in Jtor maps\")\n", "ax2.set_xlabel(r'Major radius, $R$ [m]')\n", "ax2.set_ylabel(r'Height, $Z$ [m]')\n", "\n", "eq.tokamak.plot(axis=ax2,show=False)\n", "ax2.plot(eq.tokamak.wall.R, eq.tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im2 = ax2.contourf(eq.R, eq.Z, diff, levels=np.linspace(0.01, np.max(diff), 20))\n", "cbar = plt.colorbar(im2, ax=ax2, fraction=0.09)\n", "cbar.set_label('Current [A]')\n", "\n", "# could plot difference" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the next plot, we show how to visualise the derivative of the plasma current density with respect to the metal currents (i.e. visualise the finite difference calculated Jacobian). Given we have used all the modes, the index `i` will correspond to the metal given in `eq.tokamak.coils_list[i]`.\n", "\n", "As before, just choose a mode number and plot." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extract the Jacobian (no. plasma current points x no. modes) \n", "# i.e. the full finite-difference Jacobian of plasma current density points wrt metal currents \n", "dIydI = nonlinear_solver.dIydI \n", "\n", "# choose the mode number you want to plot\n", "i = 11\n", "\n", "if i < 150:\n", " print(f\"Jacobian for metal: {eq.tokamak.coils_list[i]}\")\n", "else:\n", " print(f\"Jacobian for plasma current\")\n", "\n", "# extracts the column of the Jacobian and places into the correct grid for plotting\n", "derivs = np.full(nonlinear_solver.limiter_handler.mask_inside_limiter.shape, np.nan)\n", "derivs[nonlinear_solver.limiter_handler.mask_inside_limiter] = dIydI[:,i]\n", "derivs[derivs == 0] = np.nan # any values exactly zero should be outside LCFS\n", "\n", "# plot\n", "fig, ax = plt.subplots(1, 1, figsize=(5, 8), dpi=60)\n", "ax.grid(True, which='both')\n", "ax.set_aspect('equal')\n", "ax.set_xlim(0.1, 2.15)\n", "ax.set_ylim(-2.25, 2.25)\n", "ax.set_title(f\"Jacobian: dIy/dI[{i}]\")\n", "ax.set_xlabel(r'Major radius, $R$ [m]')\n", "ax.set_ylabel(r'Height, $Z$ [m]')\n", "plt.tight_layout()\n", "\n", "eq.tokamak.plot(axis=ax,show=False)\n", "ax.plot(eq.tokamak.wall.R, eq.tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im = ax.contourf(eq.R, eq.Z, derivs, levels=np.linspace(np.nanmin(derivs), np.nanmax(derivs), 20)) \n", "cbar = plt.colorbar(im, ax=ax, fraction=0.09)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Stability parameters\n", "\n", "We now revisit the 'Stability parameters' that were returned when we instantiated the `nl_solver` object. \n", "\n", "In addition to the deformable growth rate, we can extract other commonly used stability metrics (some that assume a deformable plasma and some that assume a rigid non-deformable plasma)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Inductive stability margin parameter\n", "\n", "The first is the inductive stability margin parameter given by equation (4) of [Humphreys et al. (2009)](https://dx.doi.org/10.1088/0029-5515/49/11/115003) and first proposed by [Portone (2005)](https://iopscience.iop.org/article/10.1088/0029-5515/45/8/021).\n", "\n", "This parameter quantifies the ability of the active coil system to stabilise vertical displacements of the plasma. It is purely inductive as it does not contain any resistance terms and is given by\n", "\n", "$$ m_s := \\lambda \\left[ -M_{m,m}^{-1} \\left( M_{m,m} + M_{m,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{I}_m} \\Bigg|_{GS} \\right) \\right], $$\n", "\n", "where the operator $\\lambda[A]$ returns the largest non-negative eigenvalue of the matrix $A$. This is the counterpart to the \"resistive\" stability margin (i.e. the deformable timescale) that is described in the \"Removing the plasma equation\" - where in this case the plasma equation has been omitted (i.e. $\\dot{I_p} = 0$).\n", "\n", "Different tokamaks have different minimum values of $m_s$ that plasmas need to remain above in order for them to be vertically controllable. Note that this parameter makes use of the Jacobian and is therefore affected by any eigenmode reduction (see below) that takes place. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(f\"Inductive stability margin = {nonlinear_solver.linearised_sol.stability_margin}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Leuer parameter\n", "\n", "The second is the stability parameter outlined in equation (6) of [Leuer (1989)](https://www.tandfonline.com/doi/abs/10.13182/FST89-A39747). See the reference for more details on its interpretation.\n", "\n", "The (generalised) Leuer stability parameter is defined as the ratio of stabilising force gradient to de-stabilising force gradient acting on the system:\n", "\n", "\\begin{equation}\n", " f := \\frac{f_{stab}}{f_{destab}} = \\frac{\\bold{I}_y^T M'_{y,m} M^{-1}_{m,m} M'_{m,y} \\bold{I}_y}{\\bold{I}_y^T M''_{y,m} \\bold{I}_m},\n", "\\end{equation}\n", "\n", "where $M'$ and $M''$ denote the first and second derivatives of the mutual inductances with respect to the $Z$ coordinate, respectively. Other matrices are defined as they are above. We can find these derivatives using the relations:\n", "\n", "\\begin{align}\n", " \\frac{d}{dZ} M_{y,m} = \\frac{d}{dZ} M (R_y, Z_y ; R_m, Z_m) &= \\frac{d}{dZ} 2 \\pi G (R_y, Z_y ; R_m, Z_m) \\\\\n", " &= \\frac{d}{dZ} 2 \\pi \\psi (R_y, Z_y ; R_m, Z_m) \\\\\n", " &= - 2 \\pi R_m B_R (R_y, Z_y ; R_m, Z_m). \n", "\\end{align}\n", "\n", "Here, $G$ is the Greens function between the plasma locations and the metal locations (which is equal to the poloidal flux $\\psi_{y,m}$ at the metal locations produced by a unit current at the plasma locations). We can then use the relation $B_R = - (d \\psi / dZ) / R$, where $B_R$ is the radial magnetic field. Note the actual calculation of this is slightly more involed when the metals are defined using filaments (see source code in nonlinear_solve.py).\n", "\n", "The second derivative is similar except we need $d B_R/ dZ$ instead of $B_R$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This stability parameter can be calculated in different ways by choosing which metals ($m$) to include in the calculation. Note, however, that no mode decomposition (see below) is performed when calculating these metrics. \n", "\n", "Also note that in this particular equilibrium the metrics look very similar because there are no currents in the passive structures!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(f\"Stabilisation force gradient (from all metals) = {nonlinear_solver.all_coils_stab_force}\")\n", "print(f\"Stabilisation force gradient (from active coils) = {nonlinear_solver.actives_stab_force}\")\n", "print(f\"Stabilisation force gradient (from passive coils) = {nonlinear_solver.passives_stab_force}\")\n", "print(\"---\")\n", "print(f\"De-stabilisation force gradient (from all metals) = {nonlinear_solver.all_coils_destab_force}\")\n", "print(f\"De-tabilisation force gradient (from active coils) = {nonlinear_solver.actives_destab_force}\")\n", "print(f\"De-tabilisation force gradient (from passive coils) = {nonlinear_solver.passives_destab_force}\")\n", "print(\"---\")\n", "print(f\"Leuer parameter (passives over actives) = {nonlinear_solver.passives_stab_force/nonlinear_solver.actives_destab_force}\")\n", "print(f\"Leuer parameter (all metals over actives) = {nonlinear_solver.all_coils_stab_force/nonlinear_solver.actives_destab_force}\")\n", "print(f\"Leuer parameter (all metals over all metals) = {nonlinear_solver.all_coils_stab_force/nonlinear_solver.all_coils_destab_force}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Back to the deformable model: Removing the plasma equation\n", "\n", "The deformable growth rates can also be computed in the case where we assume $\\bold{\\dot{I}}_p = 0$. For this, the final row/column of the state space matrices above are omitted and the system becomes:\n", "\n", "\\begin{equation*}\n", " \\left( M_{m,m} + M_{m,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{I}_m} \\Bigg|_{GS} \\right) \\dot{\\bold{I}}_m + R_{m,m} \\bold{I}_m = \\bold{V}_m - M_{m,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{\\theta}}\\Bigg|_{GS} \\dot{\\bold{\\theta}},\n", "\\end{equation*}\n", "\n", "and we find the eigenvaluesas before. \n", "\n", "This may be useful, for example, if you do not have a good intuition for what the `plasma_resistivity` parameter (in the `nl_solver` object) is. Recall that this goes into the \"lumped\" plasma equation and therefore the state space matrices when calculating the growth rates.\n", "\n", "To extract these timescales and modes, simply call the following and carry out your analysis as desired:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# accessing the growth rates (via the timescales)\n", "timescales_const_Ip = nonlinear_solver.linearised_sol.all_timescales_const_Ip # all eigenvalues: timescales (descending order) without plasma equation\n", "growth_rates_const_Ip = 1/timescales_const_Ip # all eigenvalues: growth rates (descending order) without plasma equation\n", "modes_const_Ip = nonlinear_solver.linearised_sol.all_modes_const_Ip # all eigenvectors (columns in same order as e'values) without plasma equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Mode decomposition\n", "\n", "The initial calculation of the growth rates was quite lengthy because we included all of the passive structures (each of which requires a GS solve to calculate the Jacobian matrix).\n", "\n", "In such cases where lots of passive structures are present in a tokamak, it is possible to carry out **mode decomposiiton**. This is the process by which we can exclude a number of passive structure modes that couple weakly with the plasma (thereby having a small impact on the instability timescale) and save some compute time. \n", "\n", "There are, however, a number of different ways to choose the subset of passive structure modes. These are outlined in the \"Options\" below. \n", "\n", "Before we show how to do this computationally, we outline the modifications to the original state space matrix system above in order to acheive this. Further details can be found in [Appendix 2 of Amorisco et al. (2024)](https://pubs.aip.org/aip/pop/article/31/4/042517/3286904/FreeGSNKE-A-Python-based-dynamic-free-boundary).b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### How the decomposition works\n", "\n", "To carry out the decomposition, we start by identifying the eigenvalues/vectors (characteristic frequencies and modes) of the passive structure system. We do this by diagonalising $R_{pas,pas}^{-1} M_{pas,pas}$:\n", "\n", "\\begin{equation*}\n", " R_{pas,pas}^{-1} M_{pas,pas} = P_{pas,pas} \\Lambda_{pas,pas} P_{pas,pas}^{-1},\n", "\\end{equation*}\n", "\n", "where $\\Lambda_{pas,pas}$ is the diagonal matrix of eigenvlaues and $P_{pas,pas}$ is the matrix of eigenvectors (columns). We can use this matrix to define the block diagonal change of basis matrix:\n", "\n", "\\begin{equation*}\n", " P := \\left[\\begin{array}{cc}\n", " \\mathcal{I}_{act,act} & 0 \\\\ \n", " 0 & P_{pas,pas}\n", " \\end{array}\\right],\n", "\\end{equation*}\n", "where $\\mathcal{I}_{act,act}$ is the identity matrix (with size equal to number of active coils). Supposing we retain $k$ columns of $P$ (including all of the active coils and a subset of the \"most strongly coupled\" passive structures), we can write the full vector of metal currents as \n", "\\begin{equation*}\n", " \\bold{I}_m = P \\bold{I}_d,\n", "\\end{equation*}\n", "where $\\bold{I}_d$ is the \"reduced\" set of $k$ mode currents.\n", "\n", "Now, recall the state space matrix system from above. If we apply our expression for $\\bold{I}_m$, we get:\n", "\n", "\\begin{equation*}\n", "\\begin{aligned}\n", "\\left[\\begin{array}{cc}\n", " M_{m,m}P + M_{m,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{I}_m} \\Bigg|_{GS}P & M_{m,y} \\frac{\\partial \\bold{I}_y}{\\partial I_p}\\Bigg|_{GS} \\\\ \n", "\\frac{\\hat{\\bold{I}}_y}{R_p} \\left( M_{y,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{I}_m} \\Bigg|_{GS}P + M_{y,m}P \\right) & \\frac{\\hat{\\bold{I}}_y}{R_p} M_{y,y} \\frac{\\partial \\bold{I}_y}{\\partial I_p}\\Bigg|_{GS}\n", "\\end{array}\\right]\n", "\\left[\\begin{array}{c}\n", "\\dot{\\bold{I}}_d \\\\ \n", "\\dot{I}_p\n", "\\end{array}\\right] +\n", "\n", "\\left[\\begin{array}{cc}\n", "R_{m,m}P & 0 \\\\ \n", "0 & 1\n", "\\end{array}\\right]\n", "\\left[\\begin{array}{c}\n", "\\bold{I}_d \\\\ \n", "I_p\n", "\\end{array}\\right]\n", "\n", "&=\n", "\\left[\\begin{array}{c}\n", "\\bold{V}_m \\\\ \n", "0 \n", "\\end{array}\\right]\n", "-\n", "\\left[\\begin{array}{c}\n", "M_{m,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{\\theta}}\\Bigg|_{GS} \\dot{\\bold{\\theta}} \\\\ \n", "\\frac{\\hat{\\bold{I}}_y}{R_p} M_{y,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{\\theta}}\\Bigg|_{GS} \\dot{\\bold{\\theta}}\n", "\\end{array}\\right].\n", "\n", "\\end{aligned}\n", "\\end{equation*}\n", "\n", "We end up with some modified matrices, which we can use in the same way to calculate the $k$ reduced timescales/modes. This reduces computational runtime for the Jacobian calculation as well as we are no longer calculating for all metals but a subset $k$. Note that during a linear evolutive solve (recall example notebook 5), these are the equations that will be used to evolve the mode currents (and plasma current). \n", "\n", "Next we show a few different ways to select the \"most strongly\" coupled passive structure modes. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Option 1: Fixing the number of passive structure modes explicity\n", "\n", "While the active coils will always be included as modes, you can set the number of passive structure modes directly via the `fix_n_vessel_modes` option (choose a number between 0 and the number of passive structures in the tokamak object).\n", "\n", "The process by which the passive modes are chosen is as follows:\n", "1. A first \"estimate\" of\n", " $$\n", " \\lambda_i = \\left\\| \\left( \\frac{\\partial \\mathbf{I}_y}{\\partial \\mathbf{I}_m} \\right)_{:,i} \\right\\|,\n", " $$\n", " is calculated for all passive structure modes (without fully solving the GS equation and calculating the true Jacobian). This is done by perturbing the currents in the modes (passive structures), updating the poloidal flux produced by these currents, and then evaluating the plasma current density $J_p(R,Z, \\psi)$ to obtain the perturbed $\\mathbf{I}_y$. Here, $\\psi$ is the total poloidal flux, i.e. the plasma and the metals contributions combined.\n", "\n", "2. The $\\lambda_i$ are then sorted into ascending order and with $\\hat{\\lambda} = \\max_i \\{ \\lambda_i\\}$ being the largest (strongest coupling).\n", "\n", "3. The top (strongest coupled) `fix_n_vessel_modes` are then selected for inclusion in the mode decomposition. \n", "\n", "As one might expect, increasing `fix_n_vessel_modes` will increase the timescale of the unstable mode (which in turn leads to a smaller unstable growth rate)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# instantiate the nl_solver object with fixed mode\n", "nonlinear_solver_option_1 = nonlinear_solve.nl_solver(\n", " eq=eq, \n", " profiles=profiles, \n", " GSStaticSolver=GSStaticSolver, \n", " plasma_resistivity=1e-6,\n", " fix_n_vessel_modes=50,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, we can plot how these look using similar code as before, except we need to remember that now we've done a mode decomposition, the eigenmode currents $\\bold{I}_d$ need to be transformed back into the original metal currents $\\bold{I}_m$. This is done in the plotting cell below." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# accessing the growth rates (via the timescales)\n", "timescales = nonlinear_solver_option_1.linearised_sol.all_timescales # all eigenvalues: timescales\n", "growth_rates = 1/timescales # growth rates are simply 1/timescales\n", "modes = nonlinear_solver_option_1.linearised_sol.all_modes # all eigenvectors (columns in same order as e'values)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# extracting the unstable mode\n", "mask = (timescales > 0)\n", "idx = np.where(mask)[0][0] # index of unstable mode\n", "unstable_timescales = timescales[mask]\n", "unstable_modes = np.squeeze(modes[:,mask])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# mode number (choose which one you want to visualise)\n", "# note now that there are only 63 modes instead of 151\n", "i = idx # unstable mode\n", "# i = 0 to 62 (all modes)\n", "\n", "# this function transforms the (decomposed) mode currents back to regular metal currents\n", "# (don't forget to omit the plasma current mode, i.e. the final element)\n", "mode_currents = nonlinear_solver_option_1.evol_metal_curr.IdtoIvessel(np.real(modes[0:-1,i]))\n", "\n", "# the associated instability timescale and growth rate\n", "print(f\"Mode {i} ---> {'stable' if np.real(timescales[i]) < 0 else 'unstable'}\")\n", "print(f\"Growth rate = {np.real(growth_rates[i]):.2e} [1/s]\")\n", "print(f\"Timescale = {np.real(timescales[i]):.2e} [s]\")\n", "\n", "# multiply each metal current (from the eigenvector) with its corresponding Greens matrix and sum\n", "flux = np.sum(mode_currents[:, np.newaxis, np.newaxis] * nonlinear_solver_option_1.vessel_modes_greens, axis=0)\n", "\n", "# plot\n", "fig, ax = plt.subplots(1, 1, figsize=(5, 8), dpi=60)\n", "ax.grid(True, which='both')\n", "ax.set_aspect('equal')\n", "ax.set_xlim(0.1, 2.15)\n", "ax.set_ylim(-2.25, 2.25)\n", "ax.set_title(f\"Poloidal flux (mode {i})\")\n", "ax.set_xlabel(r'Major radius, $R$ [m]')\n", "ax.set_ylabel(r'Height, $Z$ [m]')\n", "plt.tight_layout()\n", "\n", "eq.tokamak.plot(axis=ax,show=False)\n", "ax.plot(eq.tokamak.wall.R, eq.tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "im = ax.contour(eq.R, eq.Z, flux, levels=50) \n", "cbar = plt.colorbar(im, ax=ax, fraction=0.09)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Option 2: Choosing modes based on coupling and mode frequency\n", "\n", "Alternatively, one can be more precise and select modes based on how well they couple to the \"strongest\" mode (via `min_dIy_dI` and `threshold_dIy_dI`) and whether they are below some maximum mode frequency (via `max_mode_frequency`). To explain how this works, we define these parameters as follows: \n", "- $\\alpha=$ `min_dIy_dI` $\\in [0,1]$.\n", "- $\\beta=$`threshold_dIy_dI` $\\in [0,1]$ (note that $\\alpha \\leq \\beta$).\n", "- $\\omega=$`max_mode_frequency` $\\geq 0$.\n", "\n", "The approximate process to select the modes is as follows:\n", "\n", "1. Order the eigenvalues (characteristic frequencies) of \n", " $$\n", " R_{pas,pas}^{-1} M_{pas,pas},\n", " $$ \n", " in ascending order and flag any above $\\omega$ for removal. Here $R$ and $M$ are the resistance and inductance matrices for the passive structure metals only.\n", "\n", "2. As in option 1, calculate a first estimate of \n", " $$\n", " \\lambda_i = \\left\\| \\left( \\frac{\\partial \\mathbf{I}_y}{\\partial \\mathbf{I}_m} \\right)_{:,i} \\right\\|\n", " $$ \n", " for all passive structure modes without fully solving the GS equation for each (this is done for flagged modes too). This is prior to the full Jacobian calculation later on.\n", "\n", "3. Define the \"strongest\" coupled passive structure mode as $\\hat{\\lambda} = \\max_i \\{ \\lambda_i\\}$ and then keep any modes where\n", " $$\n", " \\beta \\hat{\\lambda} \\leq \\lambda_i,\n", " $$\n", " even if they had previously been flagged for removal. \n", "\n", "4. Remove any modes where\n", " $$\n", " \\lambda_i \\leq \\alpha \\hat{\\lambda},\n", " $$\n", " regardless of whether they were flagged or not.\n", "\n", "5. Any remaining modes that had been flagged and not yet removed or kept are then removed.\n", "\n", "6. For the remaining modes, the \"true\" finite-difference Jacobian is calculated and if `mode_removal=True`, any extra modes with their \"true\" $\\lambda_i \\leq \\alpha \\hat{\\lambda}$ are also removed. \n", "\n", "These selected modes are then used in the mode decomposition. \n", "\n", "The following diagram should further explain which modes are removed in the coupling-frequency space. \n", "\n", "![Mode selection diagram](data/mode_selection.jpg)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An example of how to do this follows. The print out tells you which modes are selected and which are removed at each stage of the process." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# instantiate the nl_solver object with option 2 parameters\n", "nonlinear_solver_option_2 = nonlinear_solve.nl_solver(\n", " eq=eq, \n", " profiles=profiles, \n", " GSStaticSolver=GSStaticSolver, \n", " min_dIy_dI=0.15, \n", " threshold_dIy_dI=0.5,\n", " max_mode_frequency=1e3, \n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We should note that in cases where too few modes are retained, the solver will return a message saying that the plasma is either **stable** or it is **Alfven unstable**. In the example below, we set `fix_n_vessel_modes` to zero, demonstrating what happens without any passive stabilisation. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# instantiate the nl_solver object with option 1 parameters\n", "nonlinear_solver_option_3 = nonlinear_solve.nl_solver(\n", " eq=eq, \n", " profiles=profiles, \n", " GSStaticSolver=GSStaticSolver, \n", " fix_n_vessel_modes=0, \n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Other useful optional parameters\n", "\n", "There are a number of additional optional parameters that can be selected when initialising the `nl_solver` object:\n", "- The user can choose to specify custom metal resistances (vector) and inductances (matrix) via `custom_coil_resist` and `custom_self_ind` (just ensure dimensions match the number of active and passive coils you have). \n", "- If you already calculated the Jacobian in a previous instatiation of the object, you can provide it via `dIydI` to save calculating it again. \n", "- The timestep for evolutive simulations can be set via `full_timestep` and the internal timesteps within this full one can be set with `max_internal_timestep`. If, however, `automatic_timestep` is set (as a tuple of two floats), these values are overwritten such that `full_timestep` will equal automatic_timestep[0]*unstable_mode_timescale and `max_internal_timestep` will equal automatic_timestep[1]*unstable_mode_timescale.\n", "- By setting `force_core_mask_linearization` to True, you can ensure that the finite difference calculations for the Jacobian are based on plasmas that have the same core region as each other." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### More advanced study: moving coils/passives\n", "\n", "Often of interest is to observe what happens to the growth rate parameters as you physically move an active coil or a passive structure. This can be useful for the design of new tokamaks. \n", "\n", "Here, we build a toy tokamak with three pairs of up-down symmetric active coils (P1 and P2 are connected in series, while P3 is connected in anti-series). Coils connected in anti-series are typically used for vertical control of the plasma. We will start by solving for an equilibrium and showing what it looks like. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# define the active coils \n", "active_coils = {'P1': {'1': {'R': [0.7],\n", " 'Z': [-0.8],\n", " 'dR': 0.05,\n", " 'dZ': 0.05,\n", " 'resistivity': 1.55e-08,\n", " 'polarity': 1.0,\n", " 'multiplier': 1.0},\n", " '2': {'R': [0.7],\n", " 'Z': [0.8],\n", " 'dR': 0.05,\n", " 'dZ': 0.05,\n", " 'resistivity': 1.55e-08,\n", " 'polarity': 1.0,\n", " 'multiplier': 1.0}},\n", " 'P2': {'1': {'R': [1.75],\n", " 'Z': [-0.6],\n", " 'dR': 0.05,\n", " 'dZ': 0.05,\n", " 'resistivity': 1.55e-08,\n", " 'polarity': 1.0,\n", " 'multiplier': 1.0},\n", " '2': {'R': [1.75],\n", " 'Z': [0.6],\n", " 'dR': 0.05,\n", " 'dZ': 0.05,\n", " 'resistivity': 1.55e-08,\n", " 'polarity': 1.0,\n", " 'multiplier': 1.0}},\n", " 'P3': {'1': {'R': [1.5],\n", " 'Z': [-0.5145],\n", " 'dR': 0.05,\n", " 'dZ': 0.05,\n", " 'resistivity': 1.55e-08,\n", " 'polarity': 1.0,\n", " 'multiplier': 1.0},\n", " '2': {'R': [1.5],\n", " 'Z': [0.5145],\n", " 'dR': 0.05,\n", " 'dZ': 0.05,\n", " 'resistivity': 1.55e-08,\n", " 'polarity': -1.0,\n", " 'multiplier': 1.0}}}\n", "\n", "# # define a pair of symmetric passive structures (these will have limited effect on the plasma)\n", "# passives = []\n", "# passives.append({\n", "# \"R\": [1.81,1.81,1.79,1.79],\n", "# \"Z\": [0.51,0.49,0.49,0.51],\n", "# \"name\": 'passive_lower_wall',\n", "# \"resistivity\": 5.5e-7\n", "# })\n", "# passives.append({\n", "# \"R\": [1.81,1.81,1.79,1.79],\n", "# \"Z\": [-0.51,-0.49,-0.49,-0.51],\n", "# \"name\": 'passive_upper_wall',\n", "# \"resistivity\": 5.5e-7\n", "# })\n", "\n", "# define a rectangular limiter to contain the plasma\n", "limiter = [{'R': 0.5, 'Z': -1}, {'R': 0.5, 'Z': 1}, {'R': 2, 'Z': 1}, {'R': 2, 'Z': -1}, {'R': 0.5, 'Z': -1}]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# build machine\n", "tokamak_toy = build_machine.tokamak(\n", " active_coils_data=active_coils,\n", " limiter_data=limiter,\n", " wall_data=limiter,\n", ")\n", "\n", "# build equilibrium object\n", "eq_toy = equilibrium_update.Equilibrium(\n", " tokamak=tokamak_toy,\n", " Rmin=.49, Rmax=2.01, # Radial range\n", " Zmin=-1.01, Zmax=1.01, # Vertical range\n", " nx=65, # Number of grid points in the radial direction\n", " ny=65, # Number of grid points in the vertical direction\n", ")\n", "\n", "# initialise Jtor profile object (here we use the Lao profiles) \n", "from freegsnke.jtor_update import Lao85\n", "profiles_toy = Lao85(\n", " eq=eq_toy,\n", " Ip=2e5,\n", " fvac=0.5,\n", " alpha=[58213.6],\n", " beta=[0.582136]\n", ")\n", "\n", "# assign some currents to coils \n", "eq_toy.tokamak.set_coil_current(coil_label=\"P1\", current_value=58550)\n", "eq_toy.tokamak.set_coil_current(coil_label=\"P2\", current_value=-51640)\n", "eq_toy.tokamak.set_coil_current(coil_label=\"P3\", current_value=0)\n", "\n", "# load static GS solver\n", "GSStaticSolver_toy = GSstaticsolver.NKGSsolver(eq_toy) \n", "\n", "# solve for equilibrium\n", "GSStaticSolver_toy.forward_solve(\n", " eq=eq_toy, \n", " profiles=profiles_toy, \n", " target_relative_tolerance=1e-10\n", " )\n", "\n", "# plot equilibrium\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(6, 6), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.5)\n", "eq_toy.plot(axis=ax1,show=False)\n", "tokamak_toy.plot(axis=ax1, show=False)\n", "ax1.fill(tokamak_toy.wall.R, tokamak_toy.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0)\n", "ax1.set_aspect('equal')\n", "ax1.set_xlabel(r'Major radius, $R$ $[m]$')\n", "ax1.set_ylabel(r'Height, $Z$ $[m]$')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that we have a highly-deformed plasma sitting between (but close to) the vertical stability coils. We will now see what happens to the instability timescales (and the equilbiria themselves) when we move these coils vertically away from the plasma. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# storage lists\n", "abs_P3_position = []\n", "timescales = []\n", "leuer = []\n", "eqs = []\n", "\n", "# do this for ten shifts in position\n", "for i in range(0,10):\n", "\n", " # modify P3 coil Z position\n", " active_coils['P3']['1']['Z'][0] -= 0.075 # move down 7.5cm\n", " active_coils['P3']['2']['Z'][0] += 0.075 # move up 7.5cm\n", "\n", " # build machine\n", " tokamak_toy = build_machine.tokamak(\n", " active_coils_data=active_coils,\n", " # passive_coils_data=passives,\n", " limiter_data=limiter,\n", " wall_data=limiter,\n", " )\n", "\n", " # build equilibrium object\n", " eq_toy = equilibrium_update.Equilibrium(\n", " tokamak=tokamak_toy,\n", " Rmin=.49, Rmax=2.01, # Radial range\n", " Zmin=-1.01, Zmax=1.01, # Vertical range\n", " nx=65, # Number of grid points in the radial direction\n", " ny=65, # Number of grid points in the vertical direction\n", " )\n", "\n", " # initialise Jtor profile object (here we use the Lao profiles) \n", " profiles_toy = Lao85(\n", " eq=eq_toy,\n", " Ip=2e5,\n", " fvac=0.5,\n", " alpha=[58213.6],\n", " beta=[0.582136]\n", " )\n", "\n", " # assign currents to coils \n", " eq_toy.tokamak.set_coil_current(coil_label=\"P1\", current_value=58550)\n", " eq_toy.tokamak.set_coil_current(coil_label=\"P2\", current_value=-51640)\n", " eq_toy.tokamak.set_coil_current(coil_label=\"P3\", current_value=0)\n", "\n", " # load static GS solver\n", " GSStaticSolver_toy = GSstaticsolver.NKGSsolver(eq_toy) \n", "\n", " # solve for equilibrium\n", " GSStaticSolver_toy.forward_solve(\n", " eq=eq_toy, \n", " profiles=profiles_toy, \n", " target_relative_tolerance=1e-10\n", " )\n", "\n", " # calculate the instability timescale \n", " nonlinear_solver_toy = nonlinear_solve.nl_solver(\n", " eq=eq_toy, \n", " profiles=profiles_toy,\n", " GSStaticSolver=GSStaticSolver_toy, \n", " fix_n_vessel_modes=0,\n", " )\n", "\n", " # store data\n", " abs_P3_position.append(active_coils['P3']['2']['Z'][0])\n", " timescales.append(nonlinear_solver_toy.linearised_sol.all_timescales)\n", " leuer.append(nonlinear_solver_toy.linearised_sol.all_timescales)\n", " eqs.append(eq_toy.create_auxiliary_equilibrium())\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now plot the largest timescales of the plasmas against the absolute vertical position of the P3 coil. We see that as the coils get further away from the plasma, the timescale decreases and as it passes through zero becomes Alvfen unstable (recall that the growth rate is equal to 1/timescale)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot largest timescales\n", "max_timescales = [np.max(entry) for entry in timescales]\n", "\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(6, 6), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.5)\n", "plt.plot(abs_P3_position, max_timescales, 'navy', marker='*')\n", "plt.hlines(xmin=min(abs_P3_position), xmax=max(abs_P3_position), y=0.0, linestyle='--', color='r')\n", "# ax1.set_aspect('equal')\n", "ax1.set_xlabel(r'Absolute vertical position of P3 coil $[m]$')\n", "ax1.set_ylabel(r'Largest timescale (across all modes) $[s]$')\n", "# ax1.set_yscale('log')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we plot a few of the equilibria as the coils are moved." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Indices to plot\n", "idxs = [0, 4, 9]\n", "\n", "# Create subplots\n", "fig, axes = plt.subplots(1, len(idxs), figsize=(6 * len(idxs), 6), dpi=80)\n", "\n", "# Handle single-axis case\n", "if len(idxs) == 1:\n", " axes = [axes]\n", "\n", "# Loop over equilibria and axes\n", "for ax, idx in zip(axes, idxs):\n", " plt.sca(ax)\n", " ax.grid(zorder=0, alpha=0.5)\n", " eqs[idx].plot(axis=ax, show=False)\n", " eqs[idx].tokamak.plot(axis=ax, show=False)\n", " \n", " # Fill the wall region\n", " ax.fill(\n", " eqs[idx].tokamak.wall.R,\n", " eqs[idx].tokamak.wall.Z,\n", " edgecolor='k',\n", " facecolor='w',\n", " linewidth=1.2,\n", " zorder=0\n", " )\n", " \n", " ax.set_aspect('equal')\n", " ax.set_xlabel(r'Major radius, $R$ $[m]$')\n", " ax.set_ylabel(r'Height, $Z$ $[m]$')\n", " ax.set_xlim([0.25, 2.25])\n", " ax.set_ylim([-1.5, 1.5])\n", " ax.set_title(f'Equilibrium {idx}')\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One could repeat this exercise keeping the coil positions fixed but varying a coil current. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.16" } }, "nbformat": 4, "nbformat_minor": 2 }